Non-Gaussian Velocity Distributions in Optical Lattices 
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We present a detailed experimental study of the velocity distribution of atoms cooled in an optical lattice. Our 
results are supported by full-quantum numerical simulations. Even though the Sisyphus effect, the responsible 
cooling mechanism, has been used extensively in many cold atom experiments, no detailed study of the velocity 
distribution has been reported previously. For the experimental as well as for the numerical investigation, it 
turns out that a Gaussian function is not the one that best reproduce the data for all parameters. We also fit the 
data to alternative functions, such as Lorentzians, Tsallis functions and double Gaussians. In particular, a double 
Gaussian provides a more precise fitting to our results. 

PACS numbers: 32.80.Pj, 42.50.Vk, 05.10.Ln, 05.70.Ce 



I. INTRODUCTION 

Laser cooling is now a well established technique to 
produce narrow velocity distributions for dilute samples of 
atomic gases (see e.g. Ql). The interaction between the 
atoms and the radiation modes removes kinetic energy from 
the atoms, and extremely cold samples can be obtained. In 
the standard context of Doppler or sub-Doppler laser cooling, 
atom-atom interactions are neglected and hence a thermody- 
namic temperature cannot be defined. Nevertheless, measured 
velocity distributions are generally very well fitted by a Gaus- 
sian function, and assigning a 'kinetic temperature' to the dis- 
tribution is a useful way to characterize a laser cooled atomic 
sample. 

One of the simplest theoretical models of laser cooling 
assumes a moving two-level atom interacting with counter- 
propagating pairs of laser beams, tuned slightly below the 
atomic resonance (Doppler cooling f^]). This will yield 
Doppler shifts, asymmetric with regards to velocity, and thus 
a damping force (friction). Doppler cooling is counteracted 
by momentum diffusion due to absorption and emission of 
photons. If a spatial average is taken of diffusion as well as 
friction, one obtains a stationary Gaussian velocity distribu- 
tion. This is valid since, in steady-state, most atoms have ve- 
locities well above spatial modulations in the light shift poten- 
tial (caused by the interaction between the induced dipole mo- 
ment and the light), and thus the dynamics can be described in 
terms of a Fokker-Planck equation with constant friction and 
diffusion coefficients. High irradiance results in light shifts 
of the involved energy levels that can be comparable to the 
kinetic energy, and one can no longer assume a constant ve- 
locity as atoms travel over a wavelength. Spatial averaging 
can still be performed, but one does not obtain the standard 
description of laser cooling in terms of competition between 
a friction force and a diffusion effect, since these are not sim- 
ply functions of velocity. The resulting velocity distribution 
will in this case not be Gaussian and different distributions 
have been proposed 1 3]. However, for practical Doppler cool- 
ing configurations, this effect is negligible, and there are no 
known observations of clearly non-Gaussian distributions. 



For a multilevel atom, population transfer and coherences 
between degenerate levels open up the possibility for more 
subtle cooling mechanisms. These are not limited by the ra- 
diative lifetimes of the upper levels, and can therefore lead 
to narrower distributions. In particular, Sisyphus cooling 
k4i,i^i^i7il is based on a laser beam configuration that results in 
a periodic modulation of the polarization of the light, and thus 
spatially modulated optical pumping and steady-state popula- 
tion distribution between different degenerate substates. The 
light shift will also be periodic, and will differ for different 
substates. The combination of hamiltonian motion and op- 
tical pumpingcycles transfers atomic energy to the vacuum 
modes |^ A rule of the thumb for Sisyphus cooling 

tells us that the 'temperatures' obtained correspond to kinetic 
energies that are of the order of the light shift. This behav- 
ior has been experimentally verified 19. Iiol flU [T3l down to 
kinetic temperatures of a few recoil energies. A seminal anal- 
ysis of Sisyphus cooling, by Dalibard and Cohen-Tannoudji 
[4], is again based on spatially averaged friction and diffu- 
sion coefficients. Even though the final regime corresponds 
to a situation where one can no longer assume atoms moving 
at constant velocity over many wavelengths, the scaling law 
obtained by this approach appears to be excellent. 

In mor e rig orous full quantum mechanical analyses, Castin 
et al. fis", 14] find that Sisyphus cooling ought to lead to non- 
Gaussian distributions. In particular, for irradiances close to 
the lower limit for efficient laser cooling, the effects of re- 
coils due to absorbed and emitted photons become prominent. 
Then, atomic trajectories become very irregular and the veloc- 
ity cannot be assumed to be constant. Therefore one cannot 
compute a spatially averaged velocity dependent force. Also, 
the atoms will be trapped in microscopic potential minima 
(forming optical lattices |15, 16]), and the ensemble should 
be characterized by a distribution of vibrational modes and 
unbound modes, rather than by a velocity distribution. 

Essentially all experimental investigations of Sisyphus 
cooling result in distributions that are well fitted by Gaus- 
sians. The reason for this is probably a combination of several 
facts. Many experiments are done in a regime where an aver- 
age friction coefficient seems adequate (sufficiently large light 



2 



shift). The deviations from Gaussian distributions are subtle 
and are mainly hidden in the noisy wings of the recorded dis- 
tribution. Furthermore, it is difficult to set-up an experimental 
velocity probe with the required resolution. Nevertheless, de- 
viations from Gaussian velocity distributions for laser cooled 
atoms have been reported in one recent paper iflTil . However, 
to our knowledge, there has been no systematic experimental 
study of the non-Gaussian distributions, nor any attempts to 
approach the observed distributions with more precise func- 
tions. 

In this work, we report a detailed study of velocity dis- 
tributions, as a function of the irradiance (and thus the light 
shift) for a three dimensional Sisyphus cooling configuration. 
We also perform a one-dimensional numerical simulation of 
velocity distributions, based on a full-quantum Monte-Carlo 
wave function technique. This is applied for the atomic an- 
gular momentum which is relevant in our experiment. We fit 
the recorded data, the experimental as well as the numerical, 
to different functions and compare the outcomes. 



with M being the atomic mass. In the linear regime for the 
atomic velocity, one finds 



F (v) = —av 



(2) 



In this context, a and depend on the lattice parameters 
and are independent of the velocity. d\^^ corresponds to the 

(2) 

random absorption and emission of photons while Dv rep- 
resents the fluctuations of the light-shift induced force ISj. 
The steady-state solution of Eq. Q with vanishing probability 
cuiTent {—F (v) W + AIDy (v) d^W — 0) is thus a Gaussian 
function with rms width ay = \J MDyj a: 



W (v) = Wo exp 



2MDr 



(3) 



b. Tsallis function Beyond the linear regime for atomic 
velocity, the friction force and the velocity diffusion coeffi- 
cients have to be refined into 1,13, 20.1: 



-av 



II. FITTING FUNCTIONS AND MOTIVATIONS 

The main purpose of this paper is to present more details 
about the velocity distributions of atomic samples cooled and 
trapped in optical lattices, where the Sisyphus cooling theory 
is expected to apply. A further step is to provide a function 
that gives a good approximation of the velocity distribution. 
The choice of a fitting function is made difficult by the com- 
plex dynamics of the atoms in the lattice. Indeed, even if the 
seminal process described in |4] gives very good insights in 
the dynamical behavior of the atoms, it is not sufficient in 
regimes relevant for typical experimental situations, where the 
intercombination of hamiltonian motion in the modulated po- 
tentials and optical pumping cycles, with time scales of the 
same order, makes it difficult to perform analytical calcula- 
tions flsll . Along the following lines we justify a priori the 
choice of three types of functions (Gaussian, Tsallis and dou- 
ble Gaussian) that we used to fit the experimental and the nu- 
merical recorded data. As we will see, these choices are based 
on simple considerations about well-known generalizations of 
the model presented in |4]. 

a. Gaussian function In the standard description of ID- 
Sisyphus cooling, the internal atomic state is adiabatically 
eliminated in such a way that the atomic dynamics is de- 
scribed in simple terms of a force F {v) and fluctuating forces 
of momentum diffusion coefficient Dy [v). F [v) accounts 
for the optical pumping-assisted Sisyphus cycles and Dy {v) 
corresponds, on the one hand, to the random recoils due to 
absorption and emission of photons, and on the other hand, to 
changes of potential curves. The velocity distributio n, W{v) , 
is thus governed by a Fokker-Planck equation (FPE) lll8iri9ll : 



dtW = -dy [—F{V)W]+ dy {Dy {v) ByW) ' (1) 



F{v) = - 

1 + {v/v,) 

Dy {V) = + 



(4) 



^(2) 



1 + {v/v,Y 

where v^ is the capture velocity which corresponds to the typ- 
ical atomic velocity above which the Sisyphus process breaks 
down. Now, it is straightforward to show that the steady-state 
solution with vanishing probability current of Eq. Q reads 



W{v) = Wo [I- f3{l-q)v'^ 



2MD. 



(1) 



av^ 



and (3 



a/2M 



Di'^ 



Di'^ 



(5) 



(6) 



The function in Eq. (|5} is the so-called Tsallis function and 
is in fact very general. It particularly provides a broad class 
of fitting functions including Gaussian functions (q approach- 
ing one), Lorentzian functions (q = 2) and inverted parabolas 
(q — 0). At this stage, it is interesting to note that the Tsallis 
function has been introduced in the context of non-extensive 
thermodynamics 1221 12311 . The large amount of literature in 
this context allows one to find many papers dealing with prob- 
lems already addressed in laser cooling; in particular anoma- 
lous diffusion in the presence of external forces |24, 25, "2^, 
multiplicative noise problems, and the relation to the edge of 
chaos in mixed phase space dynamics ll27ll28il . It is known 
that Sisyphus cooling can give rise to anomalous diffusion 
i29, 30], in particular for shallow optical potentials, where an 
atom can travel over many wavelengths before being trapped 
again. Even though we do not have a detailed analysis of the 
dynamics of the atoms in an optical lattice, for parameters 
corresponding to our situation, one cannot rule out anomalous 
diffusion and/or chaotic behavior 



c. Double Gaussian function As Sisyphus cooling re- 
sults in a situation where the kinetic energies of the atoms are 
of the order of the light shift potential, one can neither neglect 
atoms with lower energy ('trapped' in the potential wells) nor 
those moving more or less freely above the potential modula- 
tion (as in a 'conduction band'). This leads to a description of 
the atomic sample in terms of a bimodal dynamics. Note that 
such a bimodal description has been shown to be relevant for 
the prediction of the diffusive properties of atoms in an optical 
lattice IB li] . The kinetic equation of the 'high energy' atoms 
might very well be described by spatially averaged friction 
and diffusion coefficients resulting in a Gaussian distribution 
as shown previously. The 'low energy' atoms will be trapped, 
and subject to a different kinetic equation, and we assume that 
their velocity distribution is again a Gaussian. Our trial func- 
tion is thus the sum of two Gaussian distributions with differ- 
ent widths (double Gaussian). One corresponding to 'trapped' 
atoms and the other one to 'high energy' atoms. 



III. EXPERIMENTS 
A. Experimental setup 

The experimental setup has been described in detail previ- 
ously (see e.g. jTTl [l2il '). Briefly, we first accumulate ^■^^Cs 
atoms in a magneto-optic trap (MOT). We adjust the irradi- 
ance and the detuning, then we turn off the magnetic field and 
leave the atoms in an optical molasses with even further re- 
duced irradiance. Thus we cool the atoms to a temperature 
of 3-4 fiK. The atoms are transfered to a three-dimensional 
optical lattice, which is based on four laser beams of equal ir- 
radiance and detuning (for a review of optical lattice set-ups, 
see e.g. flsll or [16]). The detuning is a few tens of T below 
the (Fg = 4 Fe = 5) resonance for the ^^^Cs D2 line at 852 
nm (T — 2tt ■ 5.21 MHz is the linewidth of the excited state). 
The detuning (A) and irradiance (/) of the beams can be easily 
changed in order to control the depth of the light shift poten- 
tial Uq oc //|A|. The beams are aligned as in Fig. ^ two 
laser beams are linearly polarized along the a:-axis and propa- 
gate in the yz-plane symmetrically with respect to the z-axis, 
whereas the other two beams are polarized along the y-axis 
and propagate in the a;z-plane symmetrically with respect to 
z. This yields a tetragonal pattern of points with pure circular 
polarization, alternately (T+ and <t^ . These points correspond 
to potential wells where the atoms are trapped and optically 
pumped into the extreme mj? -levels (+4 and -4 respectively in 
(7^- and (T^ -wells). 

For high atomic velocities, this configuration will corre- 
spond to a three-dimensional version of the Sisyphus cooling 
model. As the atoms approach equilibrium, their kinetic en- 
ergies will get lower than the modulation depth of the optical 
potential, and thus atoms become trapped in lattice sites. They 
will get distributed in bound states, where the lowest states 
closely resemble harmonic oscillator states. 

In two different sets of runs, we let the atoms equilibrate in 
the optical lattice for 25 ms and 50 ms respectively. The veloc- 
ity distribution is then recorded with a standard time-of-flight 
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FIG. 1: Beam configuration of the 3D lin _L lin optical lattice. Two 
beam pairs propagate in the xz- and yz-planes, and are orthogonally 
polarized along the y- and x-axes respectively. They form an angle 
of 6 — 45° with the z-axis. 



method (TOF) f?*!. After the lattice period the trapping field 
is turned off, and the atoms are released in the gravitational 
field; approximately 5 cm below the trap region a thin sheet 
of resonant laser light crosses the vertical axis along which 
the atoms fall, and the induced fluorescence is recorded with 
a photo-diode. Each vertical velocity component at the time of 
release will correspond to a specific arrival time at the probe 
beam. The probe beam is carefully spatially filtered and fo- 
cussed by a cylindrical lens. The interaction region is less than 
50 /im thick, and the trapped cloud of atoms is approximately 
400 fim in diameter. This gives a velocity resolution of 0.05 
mm/s, or 0.015 ur (where iir = 3.5 mm/s is the velocity cor- 
responding to the recoil from one absorbed photon resonant 
with the D2-line). Our statistics is good enough not to con- 
tribute to this resolution. The optical lattice beams are turned 
off, by switching an acousto-optic modulator, faster than a mi- 
crosecond. This is fast enough to avoid adiabatic release of the 
atoms in the lattice, which could greatly influence the velocity 
distribution, in particular in the high velocity tails ' . 



If the optical lattice beams are turned off too slowly, the atoms may par- 
tially equilibrate in the gradually decreasing potential. There may also be 
adiabatic cooling |33]. In both these cases, the cooling during a slow turn 
off can greatly influence the velocity distribution, in particular in the high 
velocity tails. Such adiabatic switching is often used in order to achieve 
lower 'teinperatures'. 



4 



B. Experimental results 

We recorded the velocity distributions for several modula- 
tion depths and we fitted them with the functions introduced in 
section m] with a slight modification that accounts for atomic 
losses. During the long optical lattice phase, we have a con- 
stant loss of atoms, probably due to spatial diffusion. There- 
fore, the baseline is higher for atoms with a downward ve- 
locity (short times, v < 0) than it is for atoms with a up- 
ward one (v > 0). We compensate for this by adding a sharp 
step function to the fit, with the amplitude of the step as a 
free parameter. The amplitude of this step function is found 
to increase sharply for decreasing potential depths between 
Uq = 200-Er and IOOE'r. A probable reason is that spatial dif- 
fusion increases rapidly when the potential depth falls below 
some threshold, which takes place for higher potential depths 
than the threshold for cooling (usually called 'decrochage') 
j20|]. This is consistent with previous studies 113211 . In princi- 
ple, we could have used a linearly decreasing function instead 
of the step function, but then this would have had be termi- 
nated by a sharp step. We avoied this in order to minimize 
the number of free parameters and also because we wanted to 
simplify as much as possible in the absence of detailed knowl- 
edge of the loss of atoms. 

In Fig. 12] we show the rms width of the distributions, cr^, as 
a function of the depth of the optical potential Uq, as derived 
from the fits to single Gaussian functions. The width, which is 
normally associated with a kinetic temperature, increases for 
deeper potential depths as usual. 

In Figs. |3] and |4] typical recorded velocity distributions, to- 
gether with Gaussian fits, are shown for low and high mod- 
ulation depths. Figure |3l shows data taken with an equili- 
bration time of 25 ms, and for Fig. I^the equilibration time 
was 50 ms. This corresponds to typically 10^ radiative life- 
times. The plots with low irradiance are averages of twenty 
measurements and those of high irradiance of five measure- 
ments. For high values of the irradiance, a Gaussian function 
fits the velocity distribution extremely well. However, for low 
irradiance, it is clear that the wings of the distribution is not 
so well fitted. For the short equilibration time, this is more 
pronounced. 

For all data, even below 'decrochage', the attempt with 
Lorentzian fits worked very poorly. Fits to double Gaussians 
and Tsallis functions, however, reproduced recorded distribu- 
tions better than single Gaussians. In insets in figs. |3^ and 
we show fits to double Gaussians for shallow potentials. 
In Fig. |5] we compare the errors from the fits for these three 
types of functions. When the irradiance is varied, the signal- 
to-noise changes substantially, and so does the magnitude of 
the loss pedestal at short times, and the width and shape of 
the distribution. This makes it very hard to achieve a consis- 
tent normalization of the quality of the fits. The value of 
(X^ — — Xi)^, where j/i is the measured and Xi the fitted 
value) for an individual fit includes information about both 
noise and systematic deviation from the fit function, which 
are difficult to separate. The data displayed in Fig. |5]are ra- 
tios between unnormalized values of for the different fit 
functions. The displayed data are for the equilibration time of 




FIG. 2: (Color online) The rms width (av) of the measured velocity 
distributions (filled circles) as a function of the modulation depth 
of the potential. Also shown is numerically simulated data (open 
squares) in the same range (c.f. chapter lTvl . 



25 ms. The other data set has the same features. For deep po- 
tentials, all fits are essentially equally good. At more shallow 
potentials, a Tsallis function reproduces the data better than 
a Gaussian. For the whole range, a double Gaussian gives 
the best fit. For the most shallow potentials, the fitted step 
becomes too important for in order to draw any major con- 
clusion from this analysis. 

The parameter q in Eq. (|5} can be regarded as a measure of 
the shape of the distribution. A q approaching 1 will be iden- 
tical to a Gaussian distribution, whereas q — 2 corresponds to 
a Lorentzian distribution. In Fig|6] we show a plot of the fitted 
value q, for 25 ms equilibration time. For decreasing irradi- 
ances, q increases smoothly from one, and eventually reaches 
a value higher than q = 1.6. For the longer equilibration time, 
the same trend is evident, but it is much less pronounced, and 
q dose not reach higher than q = 1.3. 

The good fit to a double Gaussian can be interpreted as 
a sign of a bimodal velocity distribution. In Fig. we 
show the fitted widths of the two Gaussians for both data 
sets. This should correspond to the 'temperatures' of the two 
modes. Both these 'temperatures' increase linearly with po- 
tential depths. The areas of the two Gaussians should be a 
measure of the fraction of atoms being in one or the other 
of the modes. In Fig. 0? is the calculated relative popula- 
tions. The 'cold mode' with narrow velocity distribution al- 
ways contains most of the atoms, but the relative number of 
atoms in the 'hot mode' gets larger for decreasing potential 
depths. For potentials deeper than Uo = 250£'r there is no 
measurable portion of atoms in the 'hot modes'. The thermal 
energy of the 'hot mode' is of the same order (whithin the 
large uncertainties) as the energy barrier of the optical poten- 
tial, i.e. the modulation depth (shown in the dashed line in 
Fig. Si). 
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FIG. 3: (Color online) Experimentally recorded velocity distribu- 
tions with fits to simple Gaussians. This data is recorded with an 
equilibration time of 25 ms. For a the modulation depth of the 
optical potential was Uo = JSEr and the shown data is an aver- 
age of 20 TOF measurements. For b the corresponding facts were 
Uo = 285£'r and an average of 5 TOF measurements. The insets in 
the top right comers show magnifications of portions of the wings of 
the distributions. The inset in the top left corner of a show the same 
data with a fit to a double gaussian. 
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FIG. 4: (Color online) Experimentally recorded velocity distribu- 
tions with fits to simple Gaussians. This data is recorded with an 
equilibration time of 50 ms. For a the modulation depth of the 
optical potential was Uo = TS-Br and the shown data is an aver- 
age of 20 TOF measurements. For b the corresponding facts were 
Uo = 285i5R and an average of 5 TOF measurements. The insets in 
the top right corners show magnifications of portions of the wings of 
the distributions. The inset in the top left corner of a show the same 
data with a fit to a double gaussian. 



IV. NUMERICAL SIMULATIONS 



In order to further analyze the results of our experimen- 
tal data, we performed numerical simulations for the quantum 
dynamics of atoms undergoing Sisyphus cooling. We con- 
sider the case ofaJ = 4— >J=:5 transition, as for the ^'^^Cs 
atoms used in the experiments, and for the sake of simplicity 
we restrict the motion of the atoms into one dimension (ID). 
The laser configuration is the well-known 1D-Iin±lin config- 
uration |4] which in fact corresponds to the z-direction in our 
three dimensional (3D) experimental setup (Fig.^ with a dif- 
ferent lattice spacing. This restriction is legitimate because, 
first, the temperature has been shown to be independent of the 
lattice spacing L31t i32l] and, second, the temperature is very 



similar for both ID- and 3D- configurations (See the com- 
parison between 3D-experimental and ID-numerical results 
in Fig.|2] The slight deviation can partly be attributed to the 
difference in dimensionality). We first describe the numeri- 
cal method for the integration of the dynamics equations (sec- 
tion IIV A^ and then we present the results of the simulations 
(section llVBt . 
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FIG. 5: (Color online) Comparisons between different fits of the 
measured distribution for 25 ms equilibration time shown as ratios 
between unnormalized values of as a function of modulation 
depth of the potential. The circles are XiGauss/xIsaiiis ™d the squares 
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FIG. 6: (Color online) The fitted Tsallis g-parameter as a function of 
modulation depth of the potential for 25 ms equilibration time. 
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FIG. 7: (Color online) a) The widths of the two Gaussians as obtained 
from a fit of the data to double Gaussians, as a function of modulation 
depth of the potential. The dashed line shows the modulation depth 
in units of velocity, b) The relative population of the two modes of 
the population, obtained from the areas under the two Gaussians. For 
both a) and b), filled symbols correspond to data taken with 50 ms 
equilibration time and open symbols to 25 ms. Circles are 'temper- 
atures' and relative population of the 'hot mode', and square to the 
'cold mode'. 



A. Integration of the dynamics equations 

Consider a two level atom, with Zeeman degeneracy, inter- 
acting with a laser field 

(z, i) - ^ {e+ (z) e— * + E- [z) e+-*} . (7) 

The laser light is strong enough to be treated as a classical 
field so that we can separate the coupling between the atom 
and the electromagnetic field into a coupling to the laser field 
and a coupling to vacuum. The coupling to the laser Val in- 
duces a hamiltonian evolution for the atom. On the contrary. 



because of the coupling to the vacuum modes, Vav, the atom 
is an open quantum system for which the evolution has to be 
treated in the density matrix formalism. The evolution of the 
atom is thus governed by the generalized optical Bloch equa- 
tions (OBE) I34I I35I1 . In the regime of low saturation, where 
the experiments are performed, the excited state relaxes much 
faster than the typical evolution time of the ground state and 
thus it can be adiabatically eliminated from the OBE. The evo- 
lution of the system then reduces to a master equation for q. 
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the atomic density matrix restricted to the ground state | 



3811 : at any time. 



^ p 
where i/eff — tttt 



2M fi(A + ir/2)' 



(8) 

(9) 



Here, pis the momentum operator, A is the detuning between 
the laser and the atomic transition, AI is the mass of one atom, 

and = —D ■ aie the raising and lowering parts of 

the dipole interaction operator respectively. In Eq. (|8}, i/eff 
is a non-hermitian operator describing the atom-laser interac- 
tion^ and £,eiax is an operator describing the coupling to the 
vacuum field, i.e. spontaneous emission of photons. The in- 
tegration of the master equation is performed via a full quan- 
tum Monte-Carlo wave function method in which g 
is substituted with a set of stochastic wave functions. The 
pseudo-hamiltonian evolution (first term in Eq. (|8|i) of each 
wave function \ is governed by a Schrodinger-like equation 
involving the non-hermitian hamiltonian iJeff: 



ih- 



dt 



Hf^s I 



(10) 



Since equation ( I10> does not include the filling terms of the 
ground state from the excited state due to spontaneous emis- 
sion, \ip) is not normalized and the instantaneous spontaneous 
emission rate is given by: — ^^^^|^^y^- To take into account 
the emission of photons, the pseudo-hamiltonian evolution 
(Eq. ilO\ ) is interrupted by quantum jumps, whose repetition 
rate is determined with accordance to the spontaneous emis- 
sion rate. It follows from the emission of a photon of wave 
vector ~K and polarization ~e that the wave function is instan- 
taneously changed into 



(11) 



with relative probabilities | ."?) 



Here the excited state 



wave function \tpe) ~ 



is determined by the adia- 



?i(A+ir/2) 

batic elimination procedure of the excited state and |0) and 
111?,"?) represent the electromagnetic field states respectively 
without any photon, and with one photon of wave vector 
~K and polarization ~e. The Monte-Carlo integration then 
provides a set of time dependent stochastic wave functions 
\ip), which represent the atomic state through the average, 
a, of the density matrices associated to the wave functions, 
(7 = It is then straightforward to show that the quan- 

tum master equation for a is the same as the master equation 
for the actual density matrix g (Eq. |8}. Hence, the value of 
any observable O for the quantum system represented by g is 
equal to the ensemble average of the value of the same ob- 
servable for each stochastic wave function represented by lip) 



- The non-hermitian part of H^ff comes from the relaxation of the excited 
state. 



a 1^) = Tr (g a) . (12) 



B. Results of the simulations 

In this work, we are interested in the particular observable 
that represents the momentum distribution of the atoms. We 
have performed the full quantum Monte-Carlo integration of 
the dynamics equations for a set of 200 wave functions for var- 
ious lattice parameters (detuning and modulation depth). Be- 
cause the width of the momentum distributions are typically 
broader than several hk, the spontaneous emission pattern can 
be approximated by photons emitted along the 3D coordinate 
axes X, y or z. With such an approximation, all operators in 
Eqs. (|9j and (II 1> couple states of the form \m,p) to states 
of the form |m',p ± hk) (where m and m' represent the in- 
ternal sub-level of the atomic ground state). It is then conve- 
nient to perform the integration in the |p) -representation. The 
state is decomposed onto the basis of the \p) states (with 
p = nhk, with n an integer positive or negative). Finally, for 
usual situations considered in this work, the typical momenta 
are smaller than 20 hk, so that we take \n\ < 100. From the 
simulations, we determined the mean kinetic energy as a func- 
tion of time. After a thermalization period, the energy reaches 
a steady-state during which the momentum distribution was 
recorded and averaged. The thermalization period was chosen 
to be 1 /r, corresponding to a time in the order of a millisec- 
ond. Since the calculation is performed in ID, this time can- 
not be directly compared to the thermalization times in the 3D 
experiment. 

In order to identify whether the momentum distribution 
is compatible with a Gaussian curve or not, we first com- 
pare the root-mean-square (rms) momentum prms defined by 
Ek = P?ms/2^^ (where i?K is the mean kinetic energy of the 
atomic sample) and which represents half the width at 
1/v^ of the stationary momentum distribution. For a Gaus- 
sian distribution, those two values are equal. 

We plot in Fig.|8l the numerical results for prms and p^ as a 
function of the potential depth Uq for three different detunings 
A. We find that these values are independent of the detuning 
within the numerical errors. Several points for lower values of 
Uq have also been calculated but the atomic cloud was found 
not to thermalize. For those cases, the temperature increases 
more or less linearly and the velocity distribution becomes 
almost flat. It is also clear in Fig. |8l that p^ms and pe have 
different behaviors, pi^s reproduces the well known depen- 
dence of the kinetic energy versus the modulation depth: pi^s 
scales as \/Uq for high values of Uq and abruptly increases as 
Uq reaches very low values, typically lower than 150 (the 
point of decrochage). The minimum value of pims is found 
to be of the order of (Prms)^;^ — 4:.lhk. On the contrary, we 
find that pe increases monotonically versus Uq for low values 
as well as for high values of Uq. The minimum value of pe 
is obtained for the minimum value of C/q for which a steady- 
state velocity distribution can be obtained (Uq > 78£'r) and 
is found to be of the order of (pe)n,in — 3.4/ifc. We iden- 
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FIG. 8: Comparison between the rms momentum and the width at 
1 / \/e of the momentum distribution as a function of the potential 
depth Uo for three different detunings A = -lOF, -20r, -30r. 



tify two different regimes that can be distinguished: For Uq 
above decrochage (C/q ^ 150£'r), both pe and pmis increase 
and pe is sHghtly larger than p^-ns that is to say that the momen- 
tum distribution is wider than a Gaussian distribution with the 
same p^m?.- For Uq below decrochage ([/q ISOSr), p^ de- 
creases while prms increases rapidly as Uq decreases; the mo- 
mentum distribution has large wings and becomes narrower 
than a Gaussian distribution. These different characteristics 
are illustrated in Fig. |9] where we plot the simulated velocity 
distributions together with Gaussian fits in the two regimes 
Uo < 150Er and Uq > 150Er. 

One should note that this result is in disagreement with 
earlier calculations performed for atoms with a theoreti- 
cal J = 1/2 ^ J = 3/2 transition for which Castin et al. 
find that prms > Pe for any value of the potential depth 
Uo In fact, when running our simulation for the 

J = 1/2 ^ J — 3/2 transition, we were able to reproduce 
the results of iTsIl and we thus conclude that the discrepancy is 
due to the different atomic transitions considered in 1 13] and 
in the present work. We finally conclude that in general, the 
momentum distribution significantly differs from a Gaussian 
distribution. Moreover, we find that the threshold for p^^s at 
low values of Uq do not affect p^. 

We now turn to a more detailed analysis of the momentum 
distributions. We first fit the velocity distributions to Tsallis 
functions. The dependence of the Tsallis parameter q on the 
modulation depth is also shown in Fig. 1101 and show a linear 
dependence of q versus Uq. 

For all numerical data q differs from 1 only by less than 5 % 
and is less than 1. Moreover, q is found to tend to 1 for shallow 
potentials indicating that the best Tsallis fit is close to a Gaus- 
sian curve in disagreement with the previous discussion. The 
discrepancy between numerical simulations and experimental 
measurements may be caused by the different dimensionality 
considered in the experiments and in the simulations. 

Consider now fits to double Gaussians. We plot in Fig. II II 
numerically recorded velocity distributions in logarithmic 
scale for potential depths in the two regimes corresponding 
to shallow and deep potentials, together with fits to double 
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FIG. 9: (Color online) Numerically recorded velocity distributions 
with fits to simple Gaussians. For a, the modulation depth of the op- 
tical potential was Uo = 78-Br. For b, it was Uo ~ 235Er. The 
insets show magnifications of portions of the wings of the distribu- 
tions. 



Gaussians. For the deep potential (C/q = 293ii'R), the profile is 
essentially parabolic and thus corresponds to a Gaussian dis- 
tribution. For the shallow potential {Uq = 78Er), we clearly 
identify two contributions: in addition to a narrow parabolic 
profile (corresponding to low energetic atoms), a broad one 
(corresponding to high energetic atoms) appears. 

This supports the interpretation of the dynamics in terms 
of a bimodal atomic distribution, with each mode correspond- 
ing to 'trapped' atoms and to nearly 'free' atoms. The whole 
distribution is well fitted by a double Gaussian function. We 
plot in Fig.ll2li the widths of both the modes as functions of 
Uq. The numerical results are in good agreement with experi- 
mental ones (see Fig. 0. For shallow potentials, we find two 
Gaussian components with widths that both increase with the 
potential depth Uq, whereas for deep potentials, the 'hot com- 
ponent' is almost undetectable. Thus the route to 'decrochage' 
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FIG. 10: The g-parameter as a function of modulation depth of 
the potential, obtained from fitting the numerically computed data 
to Tsallis functions. 



for shallower potentials can be interpreted as a transfer from 
the cold mode to the hot mode. This is supported by the re- 
sults for the populations of the two Gaussian contributions to 
the velocity distribution plotted in Fig. 112b. We actually find 
that the cold mode is largely dominant even for very shallow 
potentials close to 'decrochage'. 

Finally, we compare the numerical and experimental re- 
sults. A direct quantitative comparison is not adequate, since 
the simulations are done in ID. However, qualitatively, the 
experimental data are reproduced excellently. Figures I9I12I 
show numerical data corresponding to the experimental ones 
in Figs. I3I7I The single Gaussian works for high irradiance 
but fails to fit the wings of the distribution for low irradiances. 
A Tsallis function does not fit the distribution any better than a 
single Gaussian for the numerical data. Again, the distribution 
is best fitted by a double Gaussian and this is particularly pro- 
nounced for shallow potentials. The fits to double Gaussians 
also reproduce the signature of one 'hot' and one 'cold' mode 
for shallow potentials. This strongly supports assumptions of 
a bimodal distribution. 

We find that the 'cold' mode is largely dominant even 
for very shallow potentials close to 'decrochage'in both the 
experimental (70 %) and the numerical (90 %) results (see 
Fig. 03 and 1 12b). The quantitative discrepancy between the 
limit population of the 'hot' and 'cold' modes results from the 
different dimensionnalities (3D experiments versus ID calcu- 
lations). Indeed, the fraction of bound atoms in a dD dimen- 
sion situation can be estimated by 



n 

l<fi<d 



<E„, 



(13) 



where pdoiEj) represents the population of a state of energy 
Ej , a in a normalization factor and i^max/^ stands for the max- 
imum energy of bound states in the potential wells along the 
direction /i (-Emax/^ is of the order of the potential depth). Now, 
if the space directions are separated (which is the case in the 
harmonic approximation that one can assume to be valid for 
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FIG. 1 1 : (Color online) Numerically computed velocity distribution 
for a) Uo = 78Er and b) Uo = 293Sr together with fits to Double 
Gaussians. 



bound states), then 

l<fi<d j 

Hence, assuming that all directions are equivalent, the fraction 
of atoms in non-bound states reads 



^non bound ^ J. i \^ ^ ' ^non bound J 

l<li<d 

,(1D) 



dn, 



non bound 



(15) 



Therefore, the limit populations of the 'hot' mode at 
'decrochage' are consistent in the experiments (30 % and 
d — i) and in the simulations (10 % and d = 1). 



V. CONCLUSIONS 

In this work, we have studied the velocity distributions of 
cold atomic samples obtained by Sisyphus cooling both in ex- 
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FIG. 12: a) Widths of the two Gaussians (cold mode: squares, hot 
mode: circles) as obtained from a fit of the numerical data to Double 
Gaussians, as a function of modulation depth of the potential, b) The 
relative population of the two modes of the velocity distributions, 
obtained from the areas under the two Gaussians. The cold mode 
(squares) corresponds to the large fraction whereas the hot mode (cir- 
cles) corresponds to the small fraction. 



periments with ^^'^Cs and in full quantum numerical simula- 
tions performed for the actual 4^5 transition of ^^■^Cs. We 
stressed in particular the deviation from a Gaussian distribu- 
tion. This has already been forecasted via semi-classical as 
well as quantum simulations for a simplified 1/2 3/2 tran- 
sition showing the difference of the rms velocity v^^^ and the 
velocity corresponding to half the width at 1 of the dis- 
tribution 1 13]. We recovered such a property but with a sig- 
nificantly different behavior of the ratio v^^^/v^. This shows 



that the non-Gaussian behavior of the velocity distributions is 
certainly not a trivial effect in Sisyphus cooling. 



A. Summary of our results 

Our results (experiments as well as numerical simulations) 
show that the velocity distributions are compatible with Gaus- 
sian functions for deep enough potentials (typically for Uq 
larger than a hundred recoil energies). Note that in this case, 
the atoms are trapped in the potential wells (i.e. the kinetic 
energy of the atomic cloud is significantly smaller than the po- 
tential depth). The deviation of the velocity distribution from 
Gaussian functions become more prominent for shallow light 
shift potentials. We tested several types of functions to better 
fit the shape of the velocity distributions in the range of pa- 
rameters corresponding to deeper potentials. We found that 
a better fit (coiTesponding to smaller y^) can be obtained by 
using a Tsallis function or a double Gaussian. 

Tsallis functions - The use of a Tsallis function is related 
to the details of the dynamics of atoms cooled by the Sisy- 
phus mechanism which is known to be slightly more compli- 
cate than a Brownian motion. The Tsallis function introduces 
a new parameter q which deviation from 1 measures the de- 
viation of the velocity distribution from a Gaussian function. 
The parameter q can be calculated in the 'jumping regime' 
1*39] and it is straightforward to show that q tends to 1 for high 
values of the potential depth (thus corresponding to a weak 
deviation from a Gaussian), and increase for shallow poten- 
tials. An ab initio calculation of q is more tricky in the 'os- 
cillating regime' which coiTespond to the domain of parame- 
ters for shallow potentials, near the point of 'decrochage' 113^ . 
Nevertheless we can plot the value of q corresponding to the 
best fit of the measured velocity distributions as a function of 
modulation depth. For large modulation depths, we find that 
q approaches 1, which corresponds to a Gaussian distribution, 
in agreement with the analytical calculation (see sectionHIland 
1I2II). When reducing the potential depth, we clearly observed 
an increase in q and this corresponds to a velocity distrubution 
with wings larger than in a Gaussian function. In our case the 
maximum q is close to 1 .6 and this corresponds in our exper- 
iments to a potential depth C/q — BOE'r. For Uq < 6Gi?R, the 
atomic cloud does not reach a steady state and the optical lat- 
tice disintegrates. It is interesting to note that the rms velocity 
of Tsallis distributions with q above ^ci- = 5/3 diverge |40|. If 
one would plot rms 'temperatures' of the atoms using the rms 
velocity, this would correspond to a diverging temperature. As 
one is often limited by noise in the wings of the velocity dis- 
tribution, one has a tendency to restrict the analysis to atoms 
with velocities several times below the 1/e value of the distri- 
bution. Any divergence is hence avoided. Note also that such 
divergences are very familiar: the wings of a Lorentz distri- 
bution are also known to cause a divergence of the rms value 
of the distribution. One can also recall that in the case of nar- 
row line cooling, the rms velocity diverges fTlll^?] when one 
approaches the atomic resonance, and that for very small de- 
tunings one can no longer even have a normalized distribution 
function Ii41il. 
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Double Gaussian functions - Fitting the recorded (exper- 
imental or numerical) distribution functions to double Gaus- 
sians works even better than the Tsallis function. On the one 
hand, it is not surprising that a fitting procedure with more free 
parameters gives better fits. On the other hand, the velocity 
distribution in logarithmic scale in Fig. II Ih clearly exhibits 
two components with very different widths. For deep po- 
tentials, one recovers a Gaussian distribution of 'cold atoms' 
bound in in the potential wells as expected from the above 
discussion. When decreasing the potential depth, a fraction of 
'hot atoms' grows up (for Uq < 120i?R). These atoms have an 
energy larger than the potential depth and are not trapped in 
the potential wells. We found that the fraction of 'hot atoms' 
can be significant for very shallow potentials. It reaches 30 % 
in 3D experiments and 10 % in ID numerical simulations just 
above 'decrochage' at Uq ~ 60£'r (the discrepancy between 
the experiments and the simulations is due to the different di- 
mensions as shown at the end of section HVBL This result 
strongly supports assumptions that an optical lattice has a bi- 
modal velocity distribution. A straightforward interpretation 
would be that some atoms are bound at lattice sites, whereas 
others have enough energy to move around on top of the mod- 
ulated potential. An interesting results of our work is that the 
phenomenon of 'decrochage' does not correspond to a sharp 
increase of the width of the velocity distributions correspond- 
ing to each mode but to a continuous transfer from the 'cold 
mode' to the 'hot mode'. We found that when decrochage oc- 
curs, the fraction of atoms in the 'hot mode' does not exceed 
a few tens percent. 

B. Perspectives 

The results shown in this paper stronlgy suggest that the 
simple picture for Sisyphus cooling, based on a competition 
between a diffusion and a friction (see Eq.|4}, is not adequate 
to describe the 'coldest' velocity distributions. Even though 
one has to be careful before generalizing the conclusions of 
this paper to other situations of laser cooling and/or trapping, 
the existence of two velocity modes might provide a useful 
guide to understand the dynamics and limits of laser cooling. 
One can note e.g. that for shallow potentials, one has fewer 
bound states, and the fraction of atoms in the conduction band 
gets more prominent, as shown in Figs. 171 andll2l 

These atoms will experience a friction force corresponding 
to the classical Sisyphus cooling model. The route to equilib- 
rium for the bound atoms is less clear. One hypothesis |43] 
is that bound levels are uniformly 'watered' from the conduc- 
tion band, whereas high lying levels are more likely to escape. 
Thus, the route to equilibrium is not quite a competition be- 
tween cooling and heating. A drawback with this theory is 
that it would not yield Gaussian velocity distributions. How- 
ever, this theory has the advantage that the rate of equilibra- 



tion should depend linearly on irradiance, which is consistent 
with previous experiments |44, 45]. In contrast, the standard 
Sisyphus cooling theory predicts a cooling rate independent of 
irradiance 1 4 1 . An interesting experiment would be to measure 
the velocity distribution as a function of time after a sudden 
change of the light shift potential, and see if the two popula- 
tions would evolve differently. 

Figures|3land|3seem to indicate a time dependence of the 
experimentally recorded velocity distribution. However, with 
the current data set (using only the two cooling times 25 ms 
and 50 ms), and with the current experimental uncertainties, 
we cannot draw any quantitative conclusion for the time de- 
pendence of the velocity distributrion. In future work, we will 
study the velocity distribution as a function of time. 

It would also be interesting to extend the test functions used 
in this paper to a narrow-line cooling scheme, which become 
more and more used with the laser cooling of earth-alkaline 
atoms. At this stage, one can however note, that a non normal- 
ized distribution function will have as an effect that there is no 
steady state distribution and that in this case atoms will diffuse 
to large velocities. This will appear in an experiment as a leak- 
age rate of the atoms from the optical lattice. The background 
observed in our experiment become more and more dominant 
for shallow potential wells. One might expect this to have a 
contribution from a diffusion of the atoms beyond the capture 
range of the optical lattice corresponding in practice to a non- 
normalized distribution function. A detailed analysis of the 
velocity distribution of atoms in optical lattices thus appears 
as a promising tool to study new statistical effects. 

Experiments as well as full quantum simulations (in ID 
and 3D) should allow one to get new insights in the dynamics 
of such systems. Apart from the suggestions above, future 
work could e.g. focus on the phase space dynamics of atoms 
in optical lattices and of quantum transport properties of 
ultra-cold atoms or even Bose-Einstein condensates. 
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